#include "flow_do.h"
#include "TH1D.h"
#include "TH2D.h"
#include <vector>
#include <iostream>
#include "TROOT.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "TLorentzVector.h"
#include <TStopwatch.h>


// Plot everything - Resolution, Ratio, V2, V3 // 
// Edited 2011. 10. 11  //


//******** read How many files and init number and gap

const int Init_number=50;   //500;   // starting number
const int number_gap=50;    // number gap
const int N_multi=20;       // number of output files
const int End_number=1000;
float Plot_X_range[NH][2];
float Plot_Y_range[NH][2];

// set histo range of 2nd harmonics
Plot_X_range[2][0] = Init_number - 50;  // x min
Plot_X_range[2][1] = End_number + 50;   // x max
Plot_Y_range[2][0] = 0.02;				// y min
Plot_Y_range[2][1] = 0.07; 				// y max
// set histo range of 3nd harmonics
Plot_X_range[3][0] = Init_number - 50;  // x min
Plot_X_range[3][1] = End_number + 50;   // x max
Plot_Y_range[3][0] = 0.0;				// y min
Plot_Y_range[3][1] = 0.02; 				// y max



//********* Prepare Color and Marker Setting *************
int gMarker[11] = {20,21,22,23,24,25,26,27,28,29,20};
int gColor[11] = {1,2,3,4,32,6,7,8,9,kViolet,kCyan+1};



// ******** Main Code Start ******************************
void plot(){

		TFile *filename[N_multi];


		//*************** define variables and class ***************
		// variables for each file's multiplicity
		int Multiplicity[N_multi];
		double Mul[N_multi];
		// variables for v2, v3 values
		double vn_raw[ND][NH][N_multi], vn_raw_err[ND][NH][N_multi], vn_corrected[ND][NH][N_multi];
		double vn_err[ND][NH][N_multi]={}; // how to fill value?  

		double resolution[ND][NH][N_multi];
		double resolution_err[ND][NH][N_multi];
		double resolution_2sub[ND][NH][N_multi];
		double resolution_3sub[ND][NH][N_multi];
		double resol_2sub;
		double resol_3sub[3];

		TF1 *resol_dist[NH];
		// Make TGraphErrors for vn values	
		TGraphErrors *gr_vn_raw[ND][NH];
		TGraphErrors *gr_vn_corrected[ND][NH];
		TGraphErrors *gr_vn_resolution_sub2[ND][NH];
		TGraphErrors *gr_vn_resolution_sub3[ND][NH];
		// for resolution
		TGraphErrors *gr_vn_resolution[ND][NH];
		//************************************************************

		// read each output file one by one **********************
		for(int i=0; i<N_multi; i++){
				Multiplicity[i] = Init_number+i*number_gap;
				Mul[i] = Multiplicity[i];
				filename[i] = new	TFile(Form("output/output%d.root",Multiplicity[i]),"read");
				// ******* read file *********************************************


				// ******* plot loop start here ****************************************************************
				for(int ih=2; ih<NH; ih++){ // for each harmonics
					resol_dist[ih] =  
							new TF1(Form("resol_dist%02d",ih),"sqrt(pi)/2/sqrt(2)*x*exp(-x*x/4)*(TMath::BesselI(0,x*x/4) + TMath::BesselI(1,x*x/4))",0,5);

						for(int id=1; id<ND; id++){ // for each detector
								resolution[id][ih][i] =	
									((TH1D*)filename[i]->Get(Form("hResol%02d_%02d",id,ih)))->GetMean();
								resol_2sub =	
										((TH1D*)filename[i]->Get(Form("hSub2Resol%02d_%02d",id,ih)))->GetMean();
								double x_sub = resol_dist[ih]->GetX(TMath::Sqrt(resol_2sub));
								resolution_2sub[id][ih][i] = resol_dist[ih]->Eval(TMath::Sqrt(2)*x_sub);
								
								
								
//								resolution_3sub[id][ih][i] =	
//										((TH1D*)filename[i]->Get(Form("hSub3Resol%02d_%02d_%02d",id,ih,isub)))->GetMean();


								if (resolution[id][ih][i] == 0) continue;
								resolution_err[id][ih][i] = 
										((TH1D*)filename[i]->Get(Form("hResol%02d_%02d",id,ih)))->GetMeanError();
								vn_raw[id][ih][i]=((TH1D*)filename[i]->Get(Form("hVn%02d_%02d",id,ih)))->GetMean();
								vn_raw_err[id][ih][i]= 
										((TH1D*)filename[i]->Get(Form("hVn%02d_%02d",id,ih)))->GetMeanError(); // Error of vn
								vn_corrected[id][ih][i]=vn_raw[id][ih][i]/resolution[id][ih][i]; // flow value with resolution correction
			
			// ********************** Fill graph ******************************
								gr_vn_raw[id][ih] = new TGraphErrors(N_multi, Mul, vn_raw[id][ih], 0, vn_raw_err[id][ih]); // graph for raw vn
								gr_vn_corrected[id][ih]= new TGraphErrors(N_multi, Mul, vn_corrected[id][ih], 0, vn_err[id][ih]); // graph for corrected vn 
								gr_vn_resolution[id][ih] = new TGraphErrors(N_multi, Mul, resolution[id][ih],0, resolution_err[id][ih]);  // graph for resolution
								gr_vn_resolution_sub2[id][ih] = new TGraphErrors(N_multi, Mul, resolution_2sub[id][ih],0,0);
								
						} // read from root file and calculated vn value - for each type done
				}//

		} // multiplicity loop done

		// **************** every data calculration is done ***********************
		// **************** Draw plot ********************************************

		for (ih=2; ih<NH; ih++){
				for (int id=1 ; id<ND; id++){
						TCanvas *c1 = new TCanvas();
						gStyle->SetOptStat(0); gStyle->SetOptFit(0);  gStyle->SetOptTitle(1);
						gPad->SetLogx(0); gPad->SetLogy(0);
						gPad->SetGridx(0); gPad->SetGridy(1);
						// set legend location
						leg = new TLegend(0.50 ,0.18,0.83,0.33,NULL,"brNDC");
						leg->SetTextSize(0.037);leg->SetBorderSize(0);leg->SetFillStyle(0);
						// set range of X, Y axias
						hfr = new TH2F("hfr","", 40, Plot_X_range[ih][0], Plot_X_range[ih][1], 20, Plot_Y_range[ih][0], Plot_Y_range[ih][1]);
						hfr->SetTitle(Title[id]);
						hfr->Draw();

						gr_vn_raw[id][ih]->SetLineColor(gColor[1]);
						gr_vn_raw[id][ih]->SetMarkerSize(0.6);
						gr_vn_raw[id][ih]->SetMarkerStyle(gMarker[1]);
						gr_vn_raw[id][ih]->SetMarkerColor(gColor[1]);

						gr_vn_corrected[id][ih]->SetLineColor(gColor[3]);
						gr_vn_corrected[id][ih]->SetMarkerSize(0.6);
						gr_vn_corrected[id][ih]->SetMarkerStyle(gMarker[1]);
						gr_vn_corrected[id][ih]->SetMarkerColor(gColor[3]);

						gr_vn_raw[id][ih]->Draw("same,p");
						gr_vn_corrected[id][ih]->Draw("same,p");

						leg->AddEntry(gr_vn_raw[id][ih],Form("v_{%d} raw",ih),"lp");
						leg->AddEntry(gr_vn_corrected[id][ih],Form("v_{%d}",ih),"lp");

						TLine *input_line = new TLine(Plot_X_range[ih][0],input_vn[ih], Plot_X_range[ih][1],input_vn[ih]); // input v2 line for 0.2
						input_line->SetLineWidth(2.5);
						leg->AddEntry(input_line,Form("input v_{%d}",ih),"l");
						input_line->Draw();
						leg->Draw();
				} 	// Draw and Save vn value for each type done
			}// ND loop is done
		for (ih=2; ih<NH; ih++){
				for(int id=1; id<ND; id++){   // loop for iDetector - resolution plot
						TCanvas *c2 = new TCanvas();
						gStyle->SetOptStat(0); gStyle->SetOptFit(0);  gStyle->SetOptTitle(1);
						gPad->SetLogx(0); gPad->SetLogy(0);
						gPad->SetGridx(0); gPad->SetGridy(1);
						// Set legend location
						leg = new TLegend(0.55 ,0.21,0.88,0.44,NULL,"brNDC");
						leg->SetTextSize(0.037);leg->SetBorderSize(0);leg->SetFillStyle(0);
						// Set range of X, Y axias
						hfr = new TH2F("hfr","hfr", 40, Plot_X_range[ih][0],Plot_X_range[ih][1],20, 0, 1.0); 
						hfr->SetTitle(Form("Detector = %s, harmonics = %d", Title[id].Data(), ih) );
						hfr->Draw();

						gr_vn_resolution[id][ih]->SetLineColor(gColor[1]);
						gr_vn_resolution[id][ih]->SetMarkerSize(0.5);
						gr_vn_resolution[id][ih]->SetMarkerStyle(24);
						gr_vn_resolution[id][ih]->SetMarkerColor(gColor[1]);
						gr_vn_resolution[id][ih]->Draw("same,p");
						gr_vn_resolution_sub2[id][ih]->SetMarkerSize(0.5);						
						gr_vn_resolution_sub2[id][ih]->SetMarkerStyle(24);
						gr_vn_resolution_sub2[id][ih]->SetMarkerColor(gColor[3]);
						gr_vn_resolution_sub2[id][ih]->Draw("same,p");
						leg->AddEntry(gr_vn_resolution[id][ih],"resolution","p");
						leg->AddEntry(gr_vn_resolution_sub2[id][ih],"resolution_sub2","p");
						leg->Draw();
				}  	// Draw and Save resolution plot
				// iDet loop done
		} // iHarmonics loop done  
//gr_vn_raw[1][2]->Draw();
//TH2F *href = new TH2F("a", "", 100, 400, 3600, 100, 0, 1);
//href->Draw();
//gr_vn_resolution[1][2]->Draw("same");
//gr_vn_corrected[1][2]->Draw("same");

		
}
